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Abstract 

We have used the density-functional theory to locate the freezing transitions and calculate the 
values of freezing parameters for a system of long elongated molecules which interact via the Gay- 
Berne pair potential. The pair correlation functions of isotropic phase which enter in the theory as 
input informations are found from the Percus-Yevick integral equation theory. At low temperatures 
the fluid freezes directly into the smectic A phase on increasing the density. The nematic phase 
is found to stabilize in between the isotropic and smectic A phases only at high temperatures and 
high densities. These features of the phase diagram are in good agreement with the computer 
simulation results. 

PACS numbers: 61.30 Cz, 62.20 Di, 61.30Jf 



I. INTRODUCTION 



A system consisting of anisotropic molecules is known to exhibit liquid crystalline phases 
in between the isotropic liquid and the crystalline solid. The liquid crystalline phases that 
commonly occur in a system of long elongated molecules are nematic and smectic (Sm A) 
phases jjj. In the nematic phase the full translational symmetry of the isotropic fluid phase 
(denoted as R^) is maintained but the rotational symmetry 0(3) or 50(3) (depending upon 
the presence or absence of the center of symmetry) is broken. In the simplest form of the 
axially symmetric molecules the group 5*0(3) [or 0(3)] is replaced by one of the uniaxial 
symmetry D^oh or D^o- The phase possessing the R^ADooh (denoting the semi-direct product 
of the translational group R^ and the rotational group D^oh) symmetry is known as uniaxial 
nematic phase. 

The smectic liquid crystals, in general, have a stratified structure with the long axis of 
molecules parallel to each other in layers. This situation corresponds to partial breakdown 
of translational invariance in addition to the breaking of the rotational invariance. Since a 
variety of molecular arrangements are possible within each layer, a number of smectic phases 
are possible. The simplest among them is the Sm A phase. In this phase the center of mass of 
molecules in a layer are distributed as in a two dimensional fluid but the molecular axes are 
on the average along a direction normal to the smectic layer (i.e. the director fi is normal to 
the smectic layer). The symmetry of the Sm A phase is {R^ x Z)/\Daoh where R^ corresponds 
to a two dimensional liquid structure and Z for a one-dimensional periodic structure. The 
other non-chiral and non-tilted smectic phase seen in a system of long elongated molecules is 
smectic i?/i(Sm 3^) phase. In each smectic layer of the Sm phase, the director is parallel 
to the layer normal as in the Sm A phase, and there is short-range positional but long range 
bond-orientational hexagonal orders in smectic plane. The azimuthally symmetrical x-ray 
ring of Sm A is replaced by a six- fold modulated diffused pattern jj]. The phase is thus 
characterized by a D^h point group symmetry and is uniaxial like Sm A. 

All these phases including that of the isotropic liquid and the crystalline solids are char- 
acterized by the average position and orientation of molecules and by the intermolecular 
spatial and orientational correlations. The factor responsible for the existence of these dis- 
tinguishing features is the anisotropy in both the shape of molecules and the attractive forces 
between them. The relationship between the intermolecular interactions and the relative sta- 
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bility of these phases is very intriguing and not yet fully understood. For a real system one 
faces the problem of knowing accurate intermolecular interaction as a function of intermolec- 
ular separation and orientations. This is because the mesogenic molecules are so complex 
that none of the methods used to calculate interactions between molecules can be applied 
without drastic approximations. Consequently, one is forced to use the phenomenological 
descriptions, either as a straightforward model unrelated to any particular physical system 
, or as a basis for describing by means of adjustable parameters between two molecules. 
Since our primary interest here is to relate the phases formed and their properties to the 
essential molecular factor responsible for the existence of liquid crystals and not to calculate 
the properties of any real system, the use of the phenomenological potential is justified. One 
such phenomenological model which has attracted lot of attention in computer simulations 
is one proposed by Gay and Berne [sj. 

In the Gay-Berne (GB) pair potential model, the molecules are viewed as rigid units with 
axial symmetry. Each individual molecule i is represented by a center-of-mass position ri 
and an orientation unit vector ei which is in the direction of the main symmetry axis of the 
molecule. The GB interaction energy between a pair of molecules (i, j) is given by 

xi(rij,ei,ej) = 4e(ei, ej)(/?-i2 _ ^-6^ 

where 

Here Uo is a constant defining the molecular diameter, is the distance between the center 
of mass of molecules i and j and r^j = rij/\rij\ is a unit vector along the center-center vector 
Tij = Ti — Tj. cT(ei, ej, r-jj) is the distance (for given molecular orientation) at which the 
intermolecular potential vanishes and is given by 



1-x 



1 - x^(ei.ej)2 



:i.3) 



The parameter x is a function of the ratio Xo(= o'e/cTs) which is defined in terms of the 
contact distances when the particles are end-to-end (e) and side-by-side (s), 
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Xn 



Xn 
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The orientational dependence of the potential well depth is given by a product of two 
functions 

e(ei, ej, r^) = eoe''(ei, ej)e'^(ei, ej, 1-^^) (1.5) 

where the scaling parameter eo is the well depth for the cross configuration (ei.ej = rij.ei = 
Tij.ej = 0). The first of these functions 

e(ei,ej) = [l-x'(ei.ej)r^ (1-6) 

favors the parallel alignment of the particle and so aids liquid crystal formations. The second 
function has a form analogous to o"(ei, ej, f-jj), i.e. 



i-x' 



1 - x'^(ei.ej)2 



"J/ 

The parameter x' is determined by the ratio of the well depth as 



Here k' is well-depth ratio for the side-by-side and end-to-end configuration. 

The GB model contains four parameters [xq, k', /i, u) that determine the anisotropy in the 
repulsive and attractive forces in addition to two parameters (o"o, cq) that scale the distance 
and energy, respectively. Though xq measures the anisotropy of the repulsive core, it also 
determines the difference in the depth of the attractive well between the side-by-side and the 
cross configurations. Both xq and k' play important role in stabilizing the liquid crystalline 
phases. The exact role of other two parameters fi and u are not very obvious; though they 
appear to affect the anisotropic attractive forces in a subtle way. 

The phase diagram found for the system interacting via GB potential of Eqs. (1-1)- 
(1.8) exhibits isotropic, nematic and Sm B phases j4,|Jfor = 3.0, k' = 5.0, fi = 2, and 
u = 1. An island of Sm A is, however, found to appear in the phase diagram at value 
of Xq slightly greater than 3.0 |5]. The range of Sm A extends to both higher and lower 
temperatures as Xq is increased. Also increased, the isotropic-nematic (J — A^) 

transition is seen to move to lower density (and pressure) at a given temperature. Bates and 
Luckhurst have investigated the phases and phase transition for the GB potential with 
Xq = 4.4, k' = 20.0, fi = 1 and u = 1 using the isothermal-isobaric Monte-Carlo simulations. 
At low pressure they found isotropic, Sm A and Sm B phases but not the nematic phase. 



However as the pressure is increased, the nematic phase also get stabihzed and a sequence 
of / — A^-Sm A and Sm B was found. 

In this paper we use the density functional approach to examine the different phases using 
Gay-Berne potential with the parameters chosen in the ref. 0] and compare our results with 
the results found by the simulations. The paper is organized as follows. In Sec. II, we 
describe Percus-Yevick integral equation theory for the calculation of the pair correlation 
functions of the isotropic phase. We compare our results with those found by simulations. 
In Sec. Ill the density-functional formalism has been used to locate the freezing transitions 
and freezing parameters. The paper ends with a discussion given in Sec. IV. 

II. ISOTROPIC LIQUID: PAIR CORRELATION FUNCTIONS 

In the theory of freezing to be discussed in the next section the structural informations 
of isotropic phase will be used as input data. The structural information of an isotropic 
liquid is contained in the two particle density distribution p(l,2) as the single particle 
density distribution is constant independent of position and orientation. The two-particle 
density distribution p(l,2) measures the probability of finding simultaneously a molecule 
in a volume element dridOi centered at(ri, Oi) and a second molecule in a volume element 
dr2df22 at (r2, 0,2). The pair correlation function g{l, 2) is related to 2) as 



where p(i) is the single particle density distribution. Since for the isotropic fluid p(l) = 
p(2) = Pf = {N)/V, where (A^) is the average number of molecules in the volume V, 



where r = r2 — ri. In the isotropic phase p(l, 2) depends only on distance |r2 — ri| = r, the 
orientation of molecules with respect to each other and on the direction of vector r. 

The pair distribution function g{l, 2) of the isotropic fluid is of particular interest as it 
is the lowest order microscopic quantity that contains information about the translational 
and orientational structures of the system and also has direct contact with intermolecular 
(as well as with intramolecular) interactions. For an ordered phase, on the other hand, as 
shown in the next section most of the structural informations are contained in the single 






(2.2) 
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particle distribution p(l). In the density functional theory of freezing the single particle 
distribution p(l) of an ordered phase is expressed in terms of the pair correlation function 
of the isotropic fluid (see Sec III). 

The value of g{l, 2) as a function of intermolecular separation and orientation at a given 
temperature and density is found either by computer simulation or by solving the Ornstein- 
Zernike(OZ) equation 

h{l, 2) = c(l, 2)+pf J c(l, 3)/i(2, 3)d3 (2.3) 

where d3 = dr^dfls and h{l,2) = g{l,2) — 1 and c(l,2) are, respectively, the total and 
direct pair correlation functions(DCF), using a suitable closure relation such as Percus- 
Yevick(PY) integral equation , hypernetted chain (HNC) relations. Approximations are 
introduced through these closure relations 

The Percus-Yevick closure relation is written in various equivalent form. The form adopted 
here is 

c^^(l,2) = /(l,2)[^?(l,2)-c(l,2)] (2.4) 

where /(I, 2) = exp[— /3m(1, 2)] — 1 is Mayer function , (3 = {kBT)^^ and m(1, 2) is a pair 
potential of interaction. Since for the isotropic liquid DCF is an invariant pair wise function, 
it has an expansion in body fixed (BF) frame in terms of basic set of rotational invariants, 

as 

c(ri2,fli,fl2) = E ci,i,m{n2)Yi,^{n{)Yi,^{n2) (2.5) 
where ?7T,=-m. The coefficients Q^/2m(ri2) are defined as 

ci,i,m{ri2) = J c(ri2, Oi, n2)Yi:^{n,)Y*jn2)dn,dn2 (2.6) 

Expanding all the angle dependent functions in BF frame, the OZ equation reduces to a 
set of algebraic equation in Fourier space 

hi,i,mik) = ci,i,m{k) + (-l)™^ E ci,i,m{k)hi,i,Uk) (2.7) 

where the summation is over allowed values of I3. The PY closure relation is expanded in 
spherical harmonics in the body (or space) fixed frame. The pair correlation functions are 
then found by solving these self consistently 



In our earlier work j9|, llO| we considered 30 harmonics in expansion of each orientation 
dependent function (see Eq.(2.5)); i.e. the series were truncated at the value of I indices 



equal to 6. Since the accuracy of the results depends on this number and as the anisotropy 
of the shape taken here is larger than the earlier work, we considered 54 harmonics. The 
series of each orientation dependent function was truncated at the value of / indices equal 
to 8. The numerical procedure for solving Eq(2.7) is the same as discussed in ref . jlo| . 

In Fig.l we compare the values of g(r*) = l+/iooo('"*)/4vr in BF-frame having 30 and 
54 harmonic coefficients at T*(= /c^T/eo) = 1.40 and density ri{= 7rp/(TQXo/6)= 0.44 for 
Xo = 4.4, where r* = r/ao is the reduced interparticle separation. One other projection 
of the pair-correlation is shown in Fig. 2 for the same set of parameters. It is obvious from 
these figures that even for xq = 4:A one gets good results with 30 harmonics. 

In Figs. 3 and 4 we compare the values of g{r*) and 5'22o(^*)) respectively, with those 
obtained by computer simulations 0] for rj = 0.36 and T* = 1.80. From these figures we 
find that while the PY theory gives qualitatively correct results for the pair correlation 
functions, quantitatively it underestimates the orientational correlations. The PY peak in 
g{r*) (see Fig. 3) is broad and of less hight than the one found from the simulation. Though 
peak hight in g22o{f*) compare well, the oscillation at large r* is damped more in PY results 
than in the simulation data. However, as emphasized in our earlier work the PY theory is 
reasonably accurate for the GB potential at low temperatures [lol |. 



III. FREEZING TRANSITIONS 



In the density functional approach one uses the grand thermodynamic potential defined 



as 



(3.1) 



where A is the Helmholtz free energy, /Xc the chemical potential and p(x) is a singlet dis- 
tribution function, to locate the transition. It is convenient to subtract the isotropic fiuid 

n 

thermodynamic potential from W and write it as |11|| 



AVT = W -Wf = AWi + /SW2 



(3.2) 



with 



N 



r, fi) In 



. Pf 



Ap(r, O) 



(3.3) 
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and 
AW2 



^ /dri2dl2idl22Ap(ri,Oi)c(ri2,12i,02)Ap(r2,ri2) (3.4) 



N 2pf 

Here Ap(x) = p(x) — p/, where p/ is the density of the coexisting hquid. 

The minimization of AW with respect to arbitrary variation in the ordered phase density 
subject to a constraint that corresponds to some specific feature of the ordered phase leads 
to 

^^ p(ri,»i) ^ f rfr2rf^^2c(rl2,^^l,^72;P/) X 

Pf J 

Ap(r2,122) (3.5) 

where Al is Lagrange multipher which appears in the equation because of constraint imposed 
on the minimization. 

Eq.(3.5) is solved by expanding the singlet distribution p(x) in terms of the order param- 
eters that characterize the ordered structures. One can use the Fourier series and Wigner 
rotation matrices to expand p(r, f2). Thus 

p(r, 12) = po E E Qimn{G,) exp(iG,.r)L'L(12) (3.6) 

1 Imn 

where the expansion coefficients 

2/ -I- 1 /• r 
Qlmn{G,) = jdvj dnp{vM) X 

exp(-iGq.r)L>Kl2) (3.7) 

are the order parameters, Gq the reciprocal lattice vectors, po the mean number density 
of the ordered phase and D^„(f2) the generalized spherical harmonics or Wigner rotation 
matrices. Note that for a uniaxial system consisting of cylindrically symmetric molecules 
m = n = and, therefore, one has 

p(r,0) =poEEQ^9exp(^G,.r)PKcos^) (3.8) 
I q 

and 

21 + 1 r r 

Qiq = — J dnp{v,n)exp{-iG^.r)Pi{cos9) (3.9) 

where P;(cos^) is the Legendre polynomial of degree I and 9 is the angle between the 
cylindrical axis of a molecule and the director. 
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In the present calculation we consider two orientational order parameters 

Pi = ^ = (Piicose)) (3.10) 

with 1=2 and 4, one order parameter corresponding to positional order along Z axis, 

fi = Qoo{G,) = {cos{—z)) (3.11) 

{d, being the layer spacing) and one mixed order parameter that measures the coupling 
between the positional and orientational ordering and is defined as, 

1 27r 
T = -Q2o{G,) = {cos{—z)Pi{cosd)) (3.12) 

The angular bracket in above equations indicate the ensemble average. 

The following order parameter equations are obtained by using Eqs.(3.5) and (3.9): 



Pi — dzi / sinOidOiPAcosOi) explsum] (3.13) 
2d Jo Jo 

/i = — / dzi I sin9id9iCos{ — —)exp[sum] (3-14) 
2d Jo Jo d 

T = — - / dzi I sin6id6iP2(cos6i)cos( — —) explsum] (3.15) 
2d Jo Jo d 

and change in density at the transition is found from the relation 

1 fd 



1 + Ap* = — - / dzi I sin9id9ieyiY)\sum\. (3.16) 
2d Jo Jo 



Here 



sum = Ap*C°o + 2/xcos(— )C(Jo(^i) + AC2°o(^i) + W4°o(^i) + 

d 

2Tr7 - 

2tcos{^)CI,{9,) (3.17) 

and 

CloiOi) = {^^f"pfY.i\2h + lf/\2l + lfl^C,{hLl-mmAcos9{)x 

^'^ hi 

roo 

/ ci^Li{r 12) ji{Gqri2)rl2dr 12 (3.18) 
where Cg{liLl; 000) are Clebsch- Gordon coefficients and Gg = 2T:/d. 
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In the isotropic phase all the four order parameters become zero. In the nematic phase the 
orientational order parameters P2 and P4 become nonzero but the other two parameters p, 
and T remain zero. This is because the nematic phase has no long range positional order. In 
the Sm A phase all the four order parameters are nonzero showing that the system has both 
the long range orientational and positional order along one direction. Equations (3.10)-(3.16) 
are solved self consistently using the values of direct pair correlation function harmonics 
Ci^i2i{r) evaluated in the previous section at a given value of T* and 77. This calculation is 
repeated with different values of d, the interlayer spacing. By substituting these solutions in 
Eqs.(3.2)-(3.4) we find the grand thermodynamic potential difference between ordered and 
isotropic phases; i.e. 



- ^ = -Ap* + iAp*(2 + Ap*)C'^ + \{PlCl, + PlCl) + n'Cl, + 

2j:iTCl^ + T^Cl^ (3.19) 

where 

97 I 1 

Cl^, = {2L + lf'\2L' + lfl'pfY.^\^f''C,{LL%mx 

/ CLL'i{ri2)ji{Gqri2)rl2dn2 (3.20) 

%f 

At a given temperature and density a phase with lowest grand potential is taken as the 
stable phase. Phase coexistence occurs at the values of p/ that makes —AW/N — for the 
ordered and the liquid phases. The transition from nematic to the Sm A is determined by 
comparing the values of —AW/N of these two phases at a given temperature and at different 
densities. We have done this analysis at four different temperatures, i.e. at T*=1.2, 1.4, 
1.6 and 1.8. Our results are summarized in Table 1. We see that at low temperatures, i.e. 
at 7"*= 1.2 and 1.4 the isotropic hquid freezes directly into Sm A on increasing the density. 
Nematic phase is not stable at these temperatures. However, at T*=1.6 the isotropic hquid 
on increasing the density is found first to freeze into the nematic phase at rj = 0.49 and on 
further increasing the density the nematic phase transforms into Sm A phase at 77 = 0.519. 
At high temperature, T*—1.8 the transition is found to take place between isotropic and 
nematic only for the density range considered by us. 

Using the results summarized in Table 1 we draw the phase digram in Fig. 5 in the density- 
temperature plane. The full circles represent the points given in Table 1 and lines are drawn 
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TABLE I: Values of the order parameters at the transitions for the GB potential with xq = 
4.4, fc' = 20,/x = = 1. Quantities in reduced units are d* = d/ao, Pressure P* = Pcrg/eo,/ic = 
Hc/eo, and T] = npfa^xo/Q 





Transition 


77 


d* 


Ap* 




P2 


Pi 


r 


p* 




1.20 


I-Sm A 


0.390 


3.67 


0.212 


0.64 


0.99 


0.71 


0.52 


1.06 


6.51 


1.40 


I-Sm A 


0.440 


3.90 


0.164 


0.65 


0.95 


0.66 


0.54 


1.72 


10.54 


1.60 


I-N 


0.490 


0.0 


0.087 


0.0 


0.76 


0.49 


0.0 


2.61 


15.49 




N-Sm A 


0.519 


4.00 


0.127 


0.49 


0.94 


0.67 


0.42 


3.01 


17.34 


1.80 


I-N 


0.530 


0.0 


0.072 


0.0 


0.75 


0.47 


0.0 


3.57 


20.53 



to demarcate the phase boundaries. 



IV. DISCUSSIONS 



The phase diagram shown in Fig. 5 is in good qualitative agreement with the one found 
by computer simulations . At low temperature the nematic is unstable; on increasing the 
density the fluid freezes directly into Sm A phase. The nematic phase is found to stabilize in- 
between the isotropic and Sm A phases when T* >1.5. When we compare the phase diagram 
given in Fig.5 with the one reported by Bates and Luckhurst |S] we find that while A^-Sm A 
boundary is in good agreement, the I — N boundary is shifted towards lower temperatures 
compared to that given in ref.p. As a consequence the region of existence of the nematic 
phase in the phase diagram found by us is relatively narrow. For example, at P*=3.0 we 
find T^-smAlTi-N = 0.945 whereas the value found by Bates and Luckhurst p] is 0.795. 
We also find that the values of the orientational orderparameters and the change in density 
at the transitions are higher than those given in ref 6|. 

The reason for finding the I — N phase boundary at lower temperatures compared to that 
found in simulation can be traced to the inadequacy of the PY integral equation theory dis- 
cussed in Sec. II to estimate the values of orientational correlations at higher temperatures. 

n 

Moreover, as we have shown in our earlier work jiy] that the second-order density functional 
method (discussed in Sec. Ill) has the tendency to overestimate the orientational order pa- 
rameters and the change in density. The modified weighted density approximation (MWDA) 
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version of density- functional method is expected to improve the agreement. We have 
checked it by using this theory to locate the I — N transition at T*=1.6 and 1.8. We find that 
/ - transition at T*=1.6 takes place at r/ = 0.471 with Ap* = 0.042, Pa = 0.59, P4 = 0.35 
and P*=2.38. These values are in better agreement with computer simulation than those 
given in Table 1. Similarly at T* = 1.8 the I — N transition takes place aX r] = 0.498 with 
Ap* = 0.038, P2 = 0.61, P4 = 0.35 and P*=3.10 giving better agreement compared to the 
values found from the second-order density-functional theory. The dashed line in Fig. 5 
demarcating the I — N boundary represent the results found from the MWDA version of 
the density functional method. If we take the dashed line as the phase boundary between 
isotropic and nematic phases then at P*=3.0, T]\r_smA/Ti-N = 0.890, instead of 0.945. 

In the conclusions we wish to emphasize that the freezing transitions in complex fluids 
can be predicted reasonably well with the density-functional method if the values of pair 
correlation functions in the isotropic phase are accurately known. 
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FIG. 1: Pair- Correlation functions of the centre of mass g(r*) for GB fluid with parameters xq = 
4.4, A;' = 20.0, |Lt = 1 and u = 1, at rj = 0.44 and T* = 1.40. The solid and dashed curves are, 
respectively for, 30 and 54 body-fixed harmonic coefficients. These two curves are indistinguishable 

on the scale of the figure. 
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FIG. 2: Spherical harmonic coefficient 522o('"*) in the body-fixed frame. The curves are the same 
as in Fig.l. 
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FIG. 3: Pair-Correlation function of the centre of mass g{r*) for GB fluid with parameters xq = 
4.4, k' = 20.0, /X = 1 and = 1 at ry = 0.36 and T* = 1.80. The sohd curve is our PY result and 
dashed curve is the simulation result of Bates and Luckhurst [6] . 




FIG. 4: Spherical harmonic coefficient g22o{r*) in the body-fixed frame. The curves are same as in 
Fig.3. 
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FIG. 5: Phase diagram for the GB potential with parameters xq = 4.4, k' = 20.0, = 1 and = 1. 
The soUd hnes indicate the phase boundaries obtained by using the second-order density-functional 
theory. The dashed line demarcating the isotropic and nematic phase is found from the modified 
weighted density approximation version of the density-functional theory. 
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